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EXECUTIVE  SUMMARY 


This  report  provides  the  background  to  anisoplanatic  jitter  and  how  that  relates  to  Zernike  tilt  covari¬ 
ance  and  ultimately  the  path  profile.  Space  and  Naval  Warfare  Systems  Center  Pacific’s  (SSC  Pacific) 
Optical  Channel  Characterization  in  Maritime  Atmospheres  (OCCIMA)  Python  code  is  demonstrated  with 
examples  that  match  prior  work  from  past  experiments  with  anisoplanatic  jitter.  A  new  experiment  using  a 
dual  imaging  system  instead  of  a  single  imaging  system  is  presented.  And  finally,  this  report  demonstrates 
how  to  estimate  the  C%  profile  given  jitter  measurements. 
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1.  INTRODUCTION 


In  my  previous  work  [1],  I  proposed  a  modified  differential  image  motion  monitor  (DIMM),  the  Dif¬ 
ferential  Multi-image  Motion  Monitor  (DM3),  by  using  a  stereo  imaging  apparatus.  With  this  approach, 
the  theory  for  estimating  an  averaged  path  was  possible.  Additionally,  the  image  pairs  were  aligned  for 
a  more  effective  turbulence  mitigation  algorithm.  There  exists  work  related  to  this  effort  by  White  ley  [2]. 
In  their  approach,  features  arc  also  tracked  as  in  my  method  but  the  C%  path  profile  is  estimated  instead  of 
the  average  path.  This  is  a  much  better  approach  to  atmospheric  turbulence  estimation  because  more  infor¬ 
mation  is  gained  from  the  scene. 

In  this  paper,  I  will  first  provide  the  background  in  tilt  anisoplanatism,  demonstrate  an  example  using 
the  present  version  of  Space  and  Naval  Warfare  Systems  Center  Pacific’s  Optical  Channel  Characteriza¬ 
tion  in  Maritime  Atmospheres  (OCCIMA)  Python  code,  show  how  to  model  the  DM3  and  anisoplanitic 
jitter  measurements,  and  finally  demonstrate  how  the  turbulence  strength  profile  can  be  estimated  with  any 
imaging  configuration  given  jitter  measurements. 
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2.  BACKGROUND 


2.1  INDEX-OF-REFRACTION  STRUCTURE  CONSTANT  C2N 

There  is  a  rich  set  of  literature  discussing  the  parameters  of  atmospheric  turbulence  for  which  I  refer 
the  reader  [3-7].  In  short,  the  strength  of  the  turbulence  may  be  expressed  with  an  atmospheric  index-of- 
refraction  structure  constant,  C'2.  This  constant  is  proportional  to  the  variation  of  index-of-refraction  dif¬ 
ferences  within  the  atmosphere  along  the  imaging  path. 

The  assumption  that  C2  is  homogeneous  along  an  imaging  path  is  not  a  reasonable  assumption  be¬ 
cause  of  random  wind  velocities,  heat  sources,  etc.  Instead,  a  function  of  the  path  distance,  z,  is  used  to 
characterize  C2  as  a  varying  parameter:  C2  — >  C2(z). 

Many  methods  exist  that  estimate  parameters  that  arc  path  averages  of  C2{z),  which  is  an  integration 
along  the  path.  These  methods  estimate  different  atmospheric  parameters  that  arc  related  to  a  weighted 
path  average  of  C2(z)  [2]: 
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Therefore  if  the  C2  profile  can  be  estimated,  then  the  above  turbulence  parameters  can  then  be  estimated. 


2.2  ZERNIKE  CIRCLE  POLYNOMIALS 

When  imaging  through  a  turbulent  atmosphere,  the  wavefront  arriving  at  the  sensor  (or  sensors)  is  dis¬ 
torted.  One  way  to  model  this  wavefront  distortion  is  using  the  Zernike  circle  polynomials  with  Noll’s  con¬ 
vention  [8].  Since  most  of  the  warping  is  in  the  tilt  modes  [6],  I  show  only  the  first  three  Zernike  polyno¬ 
mials  in  Table  1 .  The  wavefront  (in  circular  coordinates)  can  then  be  represented  by  the  Zernike  coeffi¬ 
cients  ap. 


4>(r,9)  =  ^2aiZi(r,0).  (5) 

i=  1 


Analogous  to  signal  decomposition,  the  wavefront  can  be  decomposed  into  Zernike  coefficients: 

fo2"fo1W(r,0)Zi(r,0)rdrd0 

/(T  Jo1  z?(r,e)rdrd6 


(6) 


The  Zernike  tilt  modes,  a,,  can  be  related  to  the  measured  angular  tilt  measurement  9  (units  of  A///)  with 
a  =  2$,  where  A  is  the  wavelength  and  D  is  the  aperture  diameter  [2]. 
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Table  1.  First  three  Zernike  polynomials. 
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Figure  1.  Illustration  of  a  general  dual  aperture  configuration  with  pairs  of  sources  from  Whiteley  [9]. 

2.3  DUAL  APERTURE  CONFIGURATION 

The  general  dual  aperture  geometry  as  portrayed  by  Whiteley  is  illustrated  in  Figure  1 .  In  this  configu¬ 
ration,  there  are  two  sources,  u  and  v,  that  are  observed  by  apertures  1  and  2,  respectively.  Then-  respective 
phases  can  be  modeled  with  Zernike  polynomials 


OO 

( PuliRpi )  —  ^  ^  &uljZi(pi ) , 
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with  phase  from  source  u  through  aperture  1  and  order  i  as  au\t  and  source  v  through  aperture  2  and  order 
j.  Whiteley  calculated  the  cross-correlation  of  the  Zernike  coefficients  as 


(auliav2j)  =  (^)5/3 
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and  D\  and  D2  arc  the  aperture  diameters  for  1  and  2,  respectively.  The  terms  nr ,  rit,  arc  Zernike  modes 
and  the  subscript  i  (or  j )  is  the  “Noll’s  order”  for  the  different  modes.  The  magnitude  of  the  aperture  sep¬ 
aration  projected  onto  layer  l  is  S[.  .//,  [■]  is  a  Bessel  function  of  order  k  of  the  first  kind.  The  angle  (pSi  is 
the  angle  between  the  coordinate  x  and  aperture  vector  s.  In  Equation  (9),  the  assumption  is  made  that  the 
turbulence  induced  phase  is  wide-sense  stationary  and  isotropic  and  that  the  von  Kami  an  spectrum  is  valid. 

The  only  interest  in  this  paper  is  for  tilt  in  the  x  and  y  directions.  Therefore  the  values  for  m  =  rrii  = 

1  for  i.  j  =  2  or  3.  Where  i  =  2  is  the  Noll  order  for  tilt  in  the  x  direction  and  i  =  3  is  the  Noll  order  for 
tilt  in  the  y  direction.  The  F,  G,  H  coefficients  simplify  to 


F  =  4, 
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Another  assumption  is  made  to  make  the  integrals  in  Equation  (9)  more  tractable.  With  xq  =  itDi^/Lq 
and  outer  scale  Lo  and  aperture  diameter  Equation  (9)  can  be  simplified  with  an  assumption  that 

D  <C  Lq  = =>  xq  ~  0, 
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the  cross-coiTelation  of  the  Zernike  coefficients  with  relation  to  C2(z)  is 
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Up  to  this  point,  I  haven’t  described  Auu  and  Av2i.  These  are  layer  scaling  factors  for  sources  u  and  v,  re¬ 
spectively.  With  the  apertures  on  the  same  x  —  y  plane,  these  scaling  factors  are  the  ratio  of  path  length  to 
the  layer  divided  by  the  total  path  length  to  the  source,  z/L ,  if  both  sources  are  L  distance  from  the  aper¬ 
tures. 

Assuming  I)  =  I)\  =  I)2  and  that  the  apertures  are  on  the  same  imaging  plane  with  respect  to  the 
turblence  layer  and  with  a  change  of  variables,  £  =  zi/L  =  Auu  =  Av2i,  si  =  s(£),  and  cr  =  2 s(^)/D, 
and  indicating  Zernike  tilt  coefficients  with  9U u  =  au \ ,  for  i  =  2,  or  3,  Equation  (19)  becomes 

(OunOvij)  =  16v/3T(8/3)  (y  )  D^3L  £  ^C2(L£)W:M(£;  <r),  (20) 

in  units  of  (A / D)2.  In  units  of  rad2,  the  covariance  is 

rL 

(euiiev2j)  =  210/sV3T(8/3)D~1^L  d£C2(£L)Witj(£;  a).  (21) 

Jo 

The  term  IT)  j  (£;  (?)  is  the  weighting  function  that  encapsulates  the  right-hand-side  of  Equation  (19)  and  is 
discussed  in  detail  in  the  next  section. 


2.4  WEIGHTING  FUNCTION 

The  closed  form  representation  of  the  weighting  function  IT)j(£;  a)  has  been  presented  by  Whiteley 
[2]  and  most  recently  Magee,  Das,  and  Welch  [10].  This  weight  function  is  dependent  only  on  the  scene 
geometry  and  not  on  the  atmosphere.  Thus  for  each  source  separation  and  image  (or  imaged  pair),  there 
exists  a  weight  function.  An  example  of  a  normalized  weight  curve  is  illustrated  in  Figure  2.  In  Figure  2, 
the  normalized  weight  is 


m(£;cr) 


Wfo*)/ 


(22) 
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0  0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8  0.9  1 

t 

Figure  2.  Example  of  normalized  path  length  (£)  vs.  normalized  weight  function  (w). 
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3.  ANISOPLANATIC  JITTER 


I  will  now  introduce  the  relationship  between  jitter  measurements,  Zernike  coefficient  covariance  and 
the  main  goal  of  this  work,  atmospheric  turbulence  C’2  profile.  It  is  best  portrayed  similar  to  Magee,  Das, 
and  Welch  [10]  where  there  arc  two  features  and  one  aperture  thus  9%  =  9un  =  9un •  Here  I  denote  the 
x-axis  as  the  Zernike  tilt  order  i  =  2  and  y-axis  as  Zernike  tilt  order  i  =  3  (9X  and  9y  respectively).  The 
first  and  second  features’  instantaneous  angular  tilts  are 

d(r)  =  (0x(r),0y{r))T,  (23) 

0{rm)  =  {9x(rm),0y(rm))T,  (24) 

respectively.  The  location  of  the  targets  arc  r  and  rm  which  arc  both  vectors  in  three-space  G  M3.  The  jitter 
measurements  which  is  the  tilt  anisoplanatism  arc  related  to  the  correlations  of  Zernike  tilt  coefficients 
with 

(<4,  ayf  =  dia§  (2($(r)^T(r))  -  2(0(r)0J  (rm)))  .  (25) 

The  relationship  between  the  tilt  covariance  matrix  components  (in  rad2)  can  be  related  to  C'2  with 

(0u(Ti)ev(Tj))  =  210/3V3T(8/3)D-1/3L  f1  d£CfaL)Wuvfa  a).  (26) 

Jo 

Jitter  can  be  predicted  with  known  C2  and  scene  geometry  by  using  Equation  (26)  into  Equation  (25).  This 
is  exactly  the  case  in  the  work  from  Magee,  Das,  and  Welch  [10]. 
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import  matplotlib 
import  numpy  as  np 
import  matplotlib . pyplot  as  pit 
import  scipy.io  as  sio 

from  OCCIMA . atmosphere . atmosgeo  import  geo,  atmos 

#  Target  separations  in  meters 

b  =  [0,  0.1,  0.2,  0.5,  1,  2,  4,  8] 

#  PI  a  tfo  rm  altitude  in  meters 

hp  =  [1550.,  1850.,  2750.,  4250.] 

#  Target  altitude  in  meters 

ht  =  [1250.] 

#  Range  to  target 

down.range  =  20000. 

#  Turbulence  m  ultiplie  r 

tm  =  2.0 

#  Diameter  of  pi  a  tfo  r  m  lens 

D  =  1  .  5 

wavelength  =  1.06e-6 

#  Number  of  turbulent  screens 

number.of .screens  =  200 

#  loading  Cn2  screens 

magee.data  = 

sio . loadmat (’OCCIMA/testdata/testdata_simulateMagee2003.mat’) 
jitters.p  =  np . zeros ((len(b) ,  len (hp) )). astype (float ) 
jitters  t  =  np . zeros ((len(b) ,  len(hp))) . astype (float ) 

#  Creates  most  of  the  atmospheric  model 

#  manipulate  the  model  in  the  next  section 

Atm  =  atmos (slant_range=down_range ,  ht=ht[0],  diameter= [D] , 
number.of  _screens  =  number_of .screens ,  wavelength=wavelength) 


Listing  1 .  OCCIMA  code  example  for  setting  up  variables  to  simulate  work  by  Magee,  Das,  and 
Welch  [10], 


4.  OCCIMA  CODE  EXAMPLE 

4.1  COMPARISON  WITH  MZA  ATMTOOLS 

Magee,  Das,  and  Welch  [10]  use  their  MZA  ATMTools  software  to  simulate  the  anisoplanatic  jitter.  In 
their  work,  they  calculate  anisoplanatic  jitter  for  different  target  separations  and  platform  altitudes.  They 
show  how  jitter  changes  when  these  scene  geometries  change.  Similarly,  I  have  OCCIMA  software  that 
achieves  the  same  results.  In  Listing  1, 1  demonstrate  how  to  setup  the  scene  geometries  using  OCCIMA 
modules.  The  target  separations  are  along  the  x-axis  which  is  parallel  with  the  x-axis  of  the  imaging  sys¬ 
tem.  Therefore,  all  vertical  (y-axis)  separation  is  zero.  In  Listing  2  each  jitter  is  measured  for  each  plat¬ 
form  altitude  and  target  separation.  Notice  that  there  arc  four  different  weight  vectors  computed  which 
is  in  Equation  (25).  As  shown  in  Figure  3,  there  is  excellent  agreement  between  the  OCCIMA  code  and 
MZA  ATMTools  for  anisoplanatic  jitter  simulation  in  Magee,  Das,  and  Welch  [10]. 


for  ii  in  range(®,  len(b)): 
for  jj  in  range (®,  len(hp)): 

#  Set  the  pi  a  tfo  r  m  altitude 

Atm. geo. hp  =  hp[jj] 

#  Make  sure  to  re  — build  the  internal  values  ( xi  ,  L,  etc.) 

Atm.buildlnternalValuesO 

#  Target  [  angular  separation  (  radians  )  , 

#  orientation  w .  r .  t  separation  vector] 

#  For  this  example  it  is  horizontal  ( x—  direction  ) 

target  =  [b [ i i ] /Atm . L ,  ©.®] 

#  A  single  lens  is  used 

platform  =  [Q.Q,  ®.©] 

#  compute  the  displacement  vector 

d  =  Atm . computePathDisplacement (platform ,  target) 
alpha  =  d  /  Atm. diameter 
zero_v  =  np . zeros.l ike (d) 

#  Compute  all  the  weights 

Wpp©  =  Atm . computeTiltCovariancePathWeight (zero_v , 

Atm. P.AXIS ,  Atm. P.AXIS) . flatten () 

Wtt®  =  Atm . computeTiltCovariancePathWeight (zero.v , 

Atm.  T_AXIS  ,  Atm.  T_AXIS)  .  flattenO 
Wppm  =  Atm . computeTiltCovariancePathWeight (alpha , 

Atm.  P_AXIS  ,  Atm.  P_AXIS)  .  flattenO 
Wttm  =  Atm . computeTiltCovariancePathWeight (alpha , 

Atm.  T_AXIS  ,  Atm.  T_AXIS)  .  flattenO 
Atm . Cn2  =  magee.data [ ’ data ’ ] [ ’ A11_CN2 ’ ] [©] [®]  [ :  ,  index] 
Atm.Cn2  =  tm* Atm . Cn2 [:]. reshape (— 1 , 1) 

#  Measure  jitter 

jitters.p [ii ,  jj],  jitters.t [ii , j j]  = 

Atm . computeAnisoplanatic J itter (Wpp© ,  Wppm,  Wtt®,  Wttm) 
index  =  index  +  1 


Listing  2.  Code  to  generate  jitter  measurements  for  each  platform  altitude  and  target  separation 
values. 
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Figure  3.  Comparison  of  OCCIMA  software  simulation  with  results  from  MZA  ATMTools 
(labeled  as  “True”). 


4.2  NEW  EXPERIMENT  WITH  DUAL  APERTURES 

I  can  do  another  experiment  that  is  not  done  in  Magee,  Das,  and  Welch  [10].  I  can  choose  instead  to 
have  dual  aperture  separations  (same  diameters)  and  a  source  pair  separated  by  2  meters;  hence  a  multi¬ 
imaging  motion  monitor.  The  experimental  setup  is  exactly  the  same  except  there  is  a  single  source  sepa¬ 
ration  pair  and  the  previous  separations  are  used  instead  for  the  platform  separations.  The  code  change  is 
very  simple  and  shown  in  Listing  3.  By  changing  the  separation  value  a  at  each  normalized  path  location 
£,  the  weight  computations  are  manipulated.  The  results  are  plotted  in  Figure  4. 

In  Figure  4,  it  is  interesting  to  note  that  there  is  an  exponential  increase  in  jitter  at  lower  altitudes  as 
the  separation  is  increased.  If  the  desire  is  to  measure  jitter,  higher  jitter  values  are  desirable  to  overcome 
quantization  limits;  therefore,  a  larger  separation  would  improve  a  jitter  measurement  system.  Unfortu¬ 
nately  there  are  limits  to  large  separations  because  the  source  must  be  imaged  in  both  apertures.  A  larger 
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separation  limits  the  minimum  distance  sources  can  be  imaged.  Also  note  that  having  an  aperture  pair  ob¬ 
serving  a  source  pair,  the  jitter  is  non-zero. 

This  experiment  also  demonstrates  that  a  new  imaging  geometry  is  simply  modeled  with  the  path  dis¬ 
placement  <7  which  is  dependent  on  the  source  and  aperture  geometries.  And  for  each  geometry  there  ex¬ 
ists  a  weight  vector. 

#  same  code  as  before 

#  A  single  target  target  pair  separated  by  2  meters 

target  =  [2.0/Atm.L,  0.0] 

#  A  dual  aperture  is  used  with  b  [  i  i  ]  separation  (meters  ) 

platform  =  [b[ii],  0.0] 

#  same  code  as  before 

Listing  3.  Excerpt  of  code  to  generate  jitter  measurements  for  each  platform  altitude  and  platform 
separation  values. 


Figure  4.  Experiment  with  jitter  measurements  from  varying  platform  dual  aperture  separations  and 
platform  height  values. 
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5.  APPROACH  TO  C2N  ESTIMATION 

To  estimate  C%  along  the  imaging  path,  the  approach  taken  by  Whiteley  [2]  and  that  I  take  as  well,  is 
to  discretize  the  path  into  slices  as  illustrated  in  Figure  5.  Assume  there  are  N  partitions  of  equal  size  along 
the  path  length  L,  we  may  model  turbulence  strength  partitions  with 


N 


C^L)  =  Y/C'ni  rect 


i— 1 


z-ti 

k 


(27) 


The  normalized  center  of  the  ith  partition  is  and  /,  is  the  normalized  partition  width.  I  assume  that  each 


slice  is  a  good  approximate  of  the  Cf,  in  that  local  region.  The  vector  c  E  KA  containing  each  compo 
nent  of  CL  is 


2  \T 


C  =  (C'ni)  ■  •  •  ,  C„N) 


Let  the  jitter  measurement  vector  m  E  with  K  observations  for  each  A:111  component  (in  rad")  be 


th 


f 7 ; 


mk  = 


uk 


=  (0u(r)0u(r))  -  (0U( r)0u(rm)). 


(28) 


(29) 


213/3v/3r(8/3)D~1/3L 

Using  Equations  (26),  (27)  and  (29)  the  jitter  measurement  component  is  related  to  with 

m  =  Pc,  (30) 

xN  with  each  row  E  M.N  having  each 


where  the  scene  geometry  is  encoded  in  the  matrix  P  E 
component  being  an  integration  of  the  weight  vectors, 


rii+U/ 2 

Phi  =  dZ  (Wuu(Z\  0)  -  WUU(L,  crk )) , 

Hi-k/  2 


(31) 


where  Wuu(Z'i  0)  is  the  weighting  along  the  path  with  cr  =  0  V  £.  The  subscripts  u  is  to  indicate  either  x  or 
y  direction  relative  to  the  separation  axis  and  k  to  indicate  which  measurement  geometry. 
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Figure  5.  Illustration  of  arbitrary  slices  along  an  imaging  path. 

We  now  have  a  compact  form  with  K  observations  with  N  turbulent  partitions  represented  in  linear 
form.  A  least-squares  estimate  of  the  vector  c  may  be  achieved  with 


c  =  P+m, 


I  /  rF  \ _ 1  T"1 

with  pseudo-inverse  P  =  (P  P)  P  . 


(32) 

(33) 
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Whiteley  refers  to  P+  as  the  “turbulence  profile  reconstructor”.  The  jitters  can  then  be  compared  for  good¬ 
ness  of  fit  with  m  =  Pc. 

As  an  example  of  estimating  c,  I  have  an  experiment  where  there  are  multiple  horizontal  feature  sepa¬ 
rations  at  the  target  plane  between  0.01  to  2  meters.  At  the  platform  plane  there  are  two  imagers  separated 
by  a  0.5  meters.  A  random  vector  c  was  generated  with  =  Cq  *  exp (att)  with  Co  =  1  x  10-14  and 
t  ~  N( 0, 1)  and  at  =  2.  The  jitter  for  this  configuration  is  in  Figure  6(a).  The  number  of  screens  chosen 
for  this  experiment  was  eight  screens.  The  true  random  vector  c  and  estimated  vector  c  using  Equation  (32) 
is  shown  in  Figure  6(b).  The  code  to  estimate  the  Crj  profile  is  deomonstrated  in  Listing  4. 


#  Measure  parallel  and  transverse  jitter  , 

#  j  i  1 1  e  r  s  _p  and  jitters  _  f  ,  respectively 

jitters  =  np . hstack (( j i tters.p ,  jitters.t)) 

#  Precompute  WppO ,  WttO  ,  Wppm ,  Wttm 
W [ ’ WppO ’ ]  =  WppO 

W [ ’ WttO ’ ]  =  WttO 
W [ ’ Wppm ’ ]  =  Wppm 
W [ ’ Wttm ’ ]  =  Wttm 

Cn2_est  =  Atm . e st imateCn2 _method_mza (W ,  jitters) 


Listing  4.  Excerpt  of  code  to  estimate  c  using  OCCIMA  python  modules. 


(a)  (b) 

Figure  6.  Jitter  and  atmospheric  profiling. 
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6.  CONCLUSION 


In  this  work  I  present  the  background  to  anisoplanatic  jitter  and  how  that  relates  to  Zernike  tilt  covari¬ 
ance  and  ultimately  the  path  profile.  I  then  demonstrate  OCCIMA  code  examples  that  match  prior 
work  from  Magee,  Das,  and  Welch  [10]  and  I  also  show  a  new  experiment  using  a  dual  imaging  system 
instead  of  a  single  imager.  Finally,  I  show  how  to  estimate  the  C%  profile  given  jitter  measurements. 

The  method  of  estimating  in  a  least-squares  approach  requires  the  matrix  P  to  be  well-conditioned. 
This  matrix  is  only  dependent  on  the  scene  geometry.  In  future  work  I  will  investigate  methods  of  using 
joint  multi-view  imaging  to  improve  the  condition  of  the  matrix  P,  and  hence  improve  accuracy  of  the  esti¬ 
mate  of  c.  Likewise,  I  will  investigate  the  eigen-decomposition  of  P  to  change  the  C*  slice  thickness  from 
uniform  to  varying  to  optimize  accuracy  of  estimation. 
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